rm(list = ls())

wd <- ".../Replication/"
setwd(wd)


# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# Table 1 ----------------------------------------

lottery <- read.csv("Data/baseline database.csv")
d <- subset(lottery, lottery$group!="Liberado por eleccion")

declare.baseline <- declare_ra(N = nrow(d), m=sum(d$lottery_winner)) 

# Make balance table
matrix <- matrix(NA, ncol = 5, nrow = 8)
colnames(matrix) <- c("Treatment","Control","Diference","Num. Obs.","RI p-value")
rownames(matrix) <- c("Male",
                      "Combatant",
                      "Apostolic name",
                      "Name length",
                      "Nr matches (5%)",
                      "Nr matches (10%)",
                      "Age (5% matches)",
                      "Age (10% matches)")

matrix[1,1] <- "All"
matrix[1,2] <- "All"
matrix[1,3] <- "Equal by design"
matrix[1,4] <- 195

matrix[2,1] <- "All"
matrix[2,2] <- "All"
matrix[2,3] <- "Equal by design"
matrix[2,4] <- 195

matrix[3,1] <- round(mean(d$apostoles[d$lottery_winner==1]), 3)
matrix[3,2] <- round(mean(d$apostoles[d$lottery_winner==0]),3)
matrix[3,3] <- round((mean(d$apostoles[d$lottery_winner==1]) - mean(d$apostoles[d$lottery_winner==0])),3)
matrix[3,4] <- 195
# conduct randomization inference
ri <- conduct_ri(
  formula = apostoles ~ lottery_winner,
  declaration = declare.baseline,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d,
  sims = 10000) 
matrix[3,5] <- as.numeric(summary(ri)[3])
rm(ri)

matrix[4,1] <- round(mean(d$namelength[d$lottery_winner==1]), 3)
matrix[4,2] <- round(mean(d$namelength[d$lottery_winner==0]),3)
matrix[4,3] <- round((mean(d$namelength[d$lottery_winner==1]) - mean(d$namelength[d$lottery_winner==0])),3)
matrix[4,4] <- 195
# conduct randomization inference
ri <- conduct_ri(
  formula = namelength ~ lottery_winner,
  declaration = declare.baseline,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d,
  sims = 10000)
matrix[4,5] <- as.numeric(summary(ri)[3])
rm(ri)

matrix[5,1] <- round(mean(d$matches5[d$lottery_winner==1]), 3)
matrix[5,2] <- round(mean(d$matches5[d$lottery_winner==0]),3)
matrix[5,3] <- round((mean(d$matches5[d$lottery_winner==1]) - mean(d$matches5[d$lottery_winner==0])),3)
matrix[5,4] <- 195
# conduct randomization inference
ri <- conduct_ri(
  formula = matches5 ~ lottery_winner,
  declaration = declare.baseline,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d,
  sims = 10000)
matrix[5,5]  <- as.numeric(summary(ri)[3])
rm(ri)

matrix[6,1] <- round(mean(d$matches10[d$lottery_winner==1]), 3)
matrix[6,2] <- round(mean(d$matches10[d$lottery_winner==0]),3)
matrix[6,3] <- round((mean(d$matches10[d$lottery_winner==1]) - mean(d$matches10[d$lottery_winner==0])),3)
matrix[6,4] <- 195
# conduct randomization inference
ri <- conduct_ri(
  formula = matches10 ~ lottery_winner,
  declaration = declare.baseline,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d,
  sims = 10000)
matrix[6,5] <- as.numeric(summary(ri)[3])
rm(ri)

matrix[7,1] <- round(mean(d$age5[d$lottery_winner==1], na.rm=TRUE), 3)
matrix[7,2] <- round(mean(d$age5[d$lottery_winner==0], na.rm=TRUE),3)
matrix[7,3] <- round((mean(d$age5[d$lottery_winner==1],na.rm=TRUE) - mean(d$age5[d$lottery_winner==0],na.rm=TRUE)),3)
d2 <- subset(d, is.na(d$age5)==FALSE)
matrix[7,4] <- nrow(d2)
# conduct randomization inference
declare.5 <- declare_ra(N = nrow(d2), m=sum(d2$lottery_winner)) 
ri <- conduct_ri(
  formula = age5 ~ lottery_winner,
  declaration = declare.5,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d2,
  sims = 10000)
matrix[7,5] <- as.numeric(summary(ri)[3])
rm(ri)

matrix[8,1] <- round(mean(d$age10[d$lottery_winner==1], na.rm=TRUE), 3)
matrix[8,2] <- round(mean(d$age10[d$lottery_winner==0], na.rm=TRUE),3)
matrix[8,3] <- round((mean(d$age10[d$lottery_winner==1],na.rm=TRUE) - mean(d$age10[d$lottery_winner==0],na.rm=TRUE)),3)
d3 <- subset(d, is.na(d$age10)==FALSE)
matrix[8,4] <- nrow(d3)
# conduct randomization inference
declare.10 <- declare_ra(N = nrow(d3), m=sum(d3$lottery_winner==1))
ri <- conduct_ri(
  formula = age10 ~ lottery_winner,
  declaration = declare.10,
  sharp_hypothesis = 0,
  assignment = "lottery_winner",
  data = d3,
  sims = 10000)
matrix[8,5] <- as.numeric(summary(ri)[3]) 
rm(ri)

print(xtable(matrix), file = "Output/Table 1.tex")
